LAMBERT1

METHODE DE RESOLUTION NUMERIQUE DU PROBLEME DE LAMBERT

	CONTENU : Mis à jour décembre 2004, revu sept 2011 I NOTATIONS ET DEFINITIONS II CLASSIFICATION DES CAS DE VOLS III FORMULATION UNIQUE DU PROBLEME POUR TOUS LES CAS ELLIPTIQUES IV EXISTENCE D'UN VOL ELLIPTIQUE V VOL HYPERBOLIQUE VI CONCLUSIONS VII RESOLUTION RAPIDE RESOLUTION FINE AVEC SPHERES D'INFLUENCE

Vers 1760 le mathématicien et mécanicien Lambert, exploitant les travaux de Newton, étudia le problème d'un trajet empruntant une conique joignant un point de départ donné P₁, à une destination d'arrivée P₂ donnée, en un temps fixé Dt, dans un champ de gravitation de constante connue m.

Le problème est donc très général et recouvre en particulier les problèmes modernes de mise au point des voyages interplanétaires, dans le champ de gravitation du Soleil, mais aussi certaines évolutions d'un satellite ou sonde, autour d'une planète.

Il m'a semblé judicieux, au moment ou tout les passionnés d'espace regardent vers Mars, de présenter une méthode de résolution accessible à tout étudiant possédant une calculatrice, pour peu qu'il ait quelques connaissances des mouvements orbitaux képlériens. Ce cours pourrait aussi intéresser les élèves des classes préparatoires, pour des TIPE.

La réponse au problème posé plus haut est apportée, après avoir été résolu géométriquement par le théorème qui suit.

THEOREME DE LAMBERT

La durée d'un voyage de P₁ à P₂, dans un champ de gravitation de constante m et de centre O, n'est fonction que des variables :

R₁+R₂ = OP₁+ OP₂

D =P₁P₂

a demi grand axe de la conique

Quelle que soit la durée de vol, il existe une trajectoire unique ellipse à moins d'un tour ou hyperbole.

NB : Ce problème, une fois résolu, permettra d'aborder l'optimisation de la date d'un tir interplanétaire, afin notamment de respecter certaines contraintes technologiques et de minimiser le coût en incréments de vitesse.

Nous allons donner dans les rubriques qui suivent les moyens de programmer la résolution numérique de ce problème qui n'a pas encore été résolu analytiquement.

I NOTATIONS-DEFINITIONS-RAPPELS :

1°) DONNEES INITIALES :

Dans un repère inertiel associé à J2000, héliocentrique écliptique( géocentrique équatorial), sont connus :

Le point de départ P₁, par ses coordonnées X₁, Y₁, Z₁

Le point d'arrivée P₂, par ses coordonnées X₂, Y₂, Z₂

La durée imposée Dt du voyage

Le voyage emprunte une trajectoire elliptique ou hyperbolique, et dans le cas elliptique la théorie est développée à moins d'un tour complet. Nous adapterons en fin de cours les calculs au cas d'un périple à plus d'un tour.

2°) REMARQUES PRELIMINAIRES :

Afin de bien fixer les idées et ne pas perdre de vue que l'objectif essentiel est le voyage interplanétaire, nous nous placerons dans ce cas. Le centre attractif O sera donc le soleil, P₁ la planète de départ (la Terre par exemple) et P₂ celle d'arrivée ( au choix, mais Mars la plus probable).

Dans les applications pratiques, le sens du mouvement sur la conique est important, bien que dans le problème de Lambert théorique il importe peu. Pour éviter un coût en vitesse prohibitif, il est clair qu'une sonde devra adopter le même sens de mouvement autour du soleil que les planètes, afin de profiter au mieux de l'effet de "fronde" lors de la sortie de la sphère d'influence de la planète de départ. On comprendra ainsi mieux les définitions de 3°).

3°) VOYAGE "COURT" ET VOYAGE "LONG" :

Durant le voyage sur la conique, la sonde allant de P1 à P2 voit son rayon vecteur "balayer" un angle d. Suivant la valeur de d par rapport à 180°, on distingue 2 cas, le voyage "court" où d < 180° et le voyage "long" avec d > 180°. Cette valeur de d ne préjuge pas du sens du mouvement. ( voir plus loin en 5° )

4°) NOTATIONS UNIVERSELLES :

Dans tous les types de voyages, elliptiques ou hyperboliques, nous introduisons les constantes suivantes, caractéristiques associées à la géométrie du voyage.

REMARQUES :

l > 0 voyage court, l < 0 voyage long.

Le lecteur vérifiera sans peine que

5°) SENS DU MOUVEMENT ORBITAL :

Sauf cas particuliers ( d = 0 ou p ) que nous ne traiterons pas, laissant le soin au lecteur de l'aborder comme cas limite, le mouvement dans le plan orbital est orienté par le vecteur unitaire classique W , défini par :

Ou encore avec d

NB 1 : Dans le cadre général du problème de Lambert, aucune condition ne sera imposée à ce vecteur. Mais en pratique pour un voyage interplanétaire, dont on sait que l'inclinaison orbitale est rarement forte, W sera "proche" de l'unitaire K orienté vers le nord écliptique. Ce sens, non obligatoire, est celui du bons sens, qui permet en général de profiter en partie ou pleinement de la vitesse de la planète au départ. On cherche donc une vitesse à l'infini raisonnable.

Donc

NB 2 : Cependant, on peut très bien aller d'un point à un autre, comme sur les figures ci-dessous. Vous pourrez vérifier que le vecteur W est pour ces 2 cas "planté" dans l'écran.

Voyage "long", sens contraire des planètes

Voyage "court", sens contraire des planètes

Il va de soi que W, portant le moment cinétique, est l'axe de mesure, dans le cas elliptique, des angles classiques que sont :

q l'anomalie vraie mesurée de 0 à 2p

j l'anomalie excentrique mesurée de 0 à 2p

M = j - esin j , mesurée de 0 à 2p

Rappelons dans le cas elliptique les relations classiques :

II CLASSIFICATION DES CAS DE VOL :

Donnons une classification des vols elliptiques, accompagnée des définitions des nouvelles variables intermédiaires y et F définies entre 0 et p.

Nous donnerons plus loin le cas hyperbolique, intéressant en théorie mais actuellement trop cher pour être pratiqué dans un vol réel, à moins qu'une propulsion nouvelle n'apparaisse, plus puissante qu'avec les technologies d'aujourd'hui.

1°) LA LOGIQUE DU CLASSEMENT :

On suppose adopté un sens de parcours, disons le plus commun, sensiblement dans le sens général des planètes autour du Soleil. Ce qui équivaut à dire que le vecteur W du repère périfocal PQW a une composante >0 sur l'axe nord-écliptique. On notera :

D	Un mouvement de descente ( mobile après l'apogée et avant le périgée, au point considéré)
M	Un mouvement de montée ( mobile ayant passé le périgée et avant l'apogée, au point considéré )
G-	Si l'angle balayé d < 180° vol dit "court"
G+	Si l'angle balayé d > 180° vol dit "long"
R₁< R₂	Le rayon vecteur de départ est < à celui d'arrivée
R₁> R₂	Le rayon vecteur de départ est > à celui d'arrivée

Les points P₁et P₂ peuvent donner lieu à 18 voyages différents, 12 en elliptique ( cas 1 à 12 ) et 6 en hyperbolique ( cas 13 à 18 ), répertoriés sur les figures qui suivent.

D &D	G- <==> R₁> R₂		Cas 7 --> R1>R2 G- Départ D & A+ Arrivée D & P-
	G+ <==> R₁< R₂		Cas4 --> R1<R2 G+ Départ D & A+ & P-, Arrivée D & P- & A+ '
M &M	G- <==> R₁< R₂		Cas 1 --> R1<R2 G- Départ M & P+ Arrivée M & A-
	G+ <==> R₁> R₂		Cas 8 --> R1>R2 G+ Départ M & P+ & A- Arrivée M & P+ & A-
D & M	G-	R₁< R₂	Cas 12 --> R1<R2 G- Départ D & P- Arrivée M & A-
		R₁> R₂	Cas 10 --> R1>R2 G- Départ D & A+ & P- Arrivée D & P- & A+
		R₁< R₂	Cas 5 --> R1<R2 G+ Départ D & P- Arrivée M & P+ & A-
		R₁> R₂	Cas 9 --> R1>R2 G+ Départ D & A+ & P- Arrivée M & P+ & A-
M & D	G-	R₁< R₂	Cas 2 --> R1<R2 G- Départ M & P+ Arrivée D & P- & A+
		R₁> R₂	Cas 6 --> R1>R2 G- Départ M & P+ Arrivée D & P- & A+
	G+	R₁< R₂	Cas 3 --> R1<R2 G+ Départ M & P+ Arrivée D & P- & A+
		R₁> R₂	Cas 11 --> R1>R2 G+ Départ M & P+ & A- Arrivée D & A+ & P-

2°) LA GEOMETRIE DES DIVERS CAS :

Le calcul du temps est donné dans chaque cas par une relation très spéciale. Quatre formules sont utilisées pour les 12 cas de calcul elliptique.

NB : Tous les cas du sens classique de 1 à 12 ont été vérifié, les données numériques sont dans le fichier source.

CAS 1

CAS 2

CAS 3

TYPE 1:

Ellipse

Ce sera pour nous le cas de base pour lequel nous établirons toutes les relations, à charge pour le lecteur ou l'étudiant d'en généraliser la portée ou de l'adapter dans les autres cas.

Pour les cas 1-2-3 on a:

R₁< R₂

j1 > j2

ce qui assure

CAS 11

Le cas 11 avec

R₁> R₂

j2> 2p-j1 & j1<p

assure aussi :

CAS 4

CAS 5

TYPE 2

Ellipse

Cas 4-5

R₁ < R₂

j1 > j2

ce qui assure

CAS 6

CAS 7

TYPE 3

Ellipse

Cas 4-5-6-7

R₁ > R₂

j₁ < j₂

les conditions assurent comme pour le type 1

CAS 8

CAS 9

CAS 10

CAS 12

TYPE 4

Ellipse

j₁ > j₂

Ce qui assure

CONCLUSION : Il ne reste que 2 configurations de calcul chacune pour 6 cas:

CAS

1-2-3-6-7-11

CAS

4-5-8-9-10-12

Pour un calcul du temps de parcours unique

III UNE FORMULATION UNIQUE POUR TOUS LES CAS ELLIPTIQUES:

1°) RELATIONS DE BASE EN FONCTION DE y et F. :

Nous laissons au lecteur le soin d'établir les relations qui suivent, vérifiées dans tous les cas par les variables y et F.

Seules sont indiquées les formules à utiliser par leur n°, formules rappelées précédemment.

En fait ces variables ne servent que d'intermédiaires et 2 variables définitives a et b sont introduites :

2°) VARIABLES a ET b. :

Deux dernières variables sont introduites, conduisant à des relations très "symétriques".

a = y + F

b = F - y

3°) RELATIONS DE BASE EN FONCTION DE a ET b. :

Le lecteur vérifiera les résultats suivants déduits des relations de III 1°), applicables à tous les cas elliptiques :

Le calcul de la durée du voyage s'exprime alors par l'équation ( I) :

(I)

IV EXISTENCE D'UN VOL ELLIPTIQUE :

Le temps de vol Dt étant fixé, une trajectoire "courte" ou "longue" n'existe que si une valeur de a existe, vérifiant l'équation (I) ci-dessus.

1°) ETUDE DE F(a) :

Une étude mathématique fine, montre que sur l'intervalle [0 2p] des valeurs possibles du paramètre a, la dérivée fe F est positive et tend vers 0 lorsque a tend vers 0. Donc F est croissante sur cet intervalle et son minimum est obtenu pour a = 0, qui pose une ambiguïté de limite qui se résout sans difficulté par un développement limité.

2°) EXISTENCE DU CAS ELLIPTIQUE :

On montre ainsi qu'un voyage elliptique ne peut exister que

car la fonction F(a) est monotone.

On peut encore traduire cette relation uniquement avec la géométrie du triangle OP₁P₂, par la condition :

V VOL HYPERBOLIQUE :

On "sent" bien que si la durée imposée du voyage est trop courte, la vitesse de tir devra augmenter et inévitablement jusqu'à dépasser la vitesse de libération en P1, il devra être fait usage d'une trajectoire hyperbolique. Ainsi, non limité en théorie par la vitesse de tir, un voyage sera toujours possible, y compris pour un temps de vol tendant vers 0.

Outre le fait que ce type de vol a peu de chance d'être pratiqué dans les années qui viennent, pour éviter de nouveaux calculs, nous incitons le lecteur à établir par analogie, la plupart des résultats de ce chapitre.

1°) RAPPELS SUR L'HYPERBOLE :

La figure ci-dessous rappelle les données principales de l'hyperbole :

Les asymptotes

Le repère périfocal PQW

Le foyer O

Les coordonnées X et Y dans P et Q

Donnons les principales relations, avec notamment le paramétrage utilisant le réel j, nommé ANOMALIE EXCENTRIQUE HYPERBOLIQUE et pouvant prendre toute valeur réelle possible positive ou négative.

Avec les relations utiles, mais non usuelles :

2°) CAS DE VOL HYPERBOLIQUE :

Hyperbole

R₁ < R₂

j₁ < j₂

d > p

Hyperbole

R₁ > R₂

j₁ < j₂

d < p

Hyperbole

R₁ > R₂

j₁ < j₂

d > p

Hyperbole

R₁ < R₂

j₁ < j₂

d < p

3°) CONVENTIONS :

Les constantes géométriques S, D, l, d sont définies de la même manière que dans le cas elliptique. Il en est de même du vecteur unitaire W orientant le mouvement orbital. La figure qui suit résume les données.

On appelle j₁ et j₂ , avec j₁ < j₂, les anomalies excentriques hyperboliques associées à P₁ et P₂.

Comme pour le cas elliptique, on introduit de nouvelles variables F et Y, puis a et b, définies par :

La constante l conserve la même définition que pour le cas elliptique

3°) FORMULES TRANSITOIRES AVEC Y et F :

Le lecteur établira, par des calculs analogues à ceux du cas elliptique, et ceci quel que soit le cas hyperbolique, les relations suivantes:

4°) FORMULES DECISIVES AVEC a et b :

Enfin, pour achever les calculs, vous montrerez qu'une seule relation suffit pour tous les cas hyperboliques de 4 à 7 :

l pouvant être naturellement positif ou négatif, suivant la valeur de d.

5°) ETUDE DE G(a) : :

Le lecteur démontrera ou admettra que la fonction G(a) est décroissante sur R+, son maximum coïncidant avec le minimum de F(a). La valeur minimale de G(a) est 0 quand a tend vers l'infini.

VI CONCLUSIONS :

Nous constatons donc qu'un voyage à temps fixé, à moins d'un tour pour le cas elliptique, est toujours possible quel que soit ce temps, soit sur une ellipse, soit sur une hyperbole. Le problème revient alors à résoudre, pour chaque cas, une équation donnant a.

Le reste du calcul c'est à dire la résolution conduisant aux valeurs des anomalies excentriques puis notamment au calcul des paramètres orbitaux, relève des autres cours déjà présentés et fera l'objet d'un excellent projet pour des étudiants de DESS.

REMARQUE FINALE :

Dans le cas elliptique a peut varier de 0 à 2p et pour le cas hyperbolique de 0 l'infini.

Une variable X appelée variable de Lambert, peut être introduite, en posant :

Cette variable X caractérise complètement la nature du voyage.

VI ADAPTATION AU VOYAGE ELLIPTIQUE A PLUS D'UN TOUR :

1°) ADAPTATION DES CALCULS :

Le lecteur adaptera les calculs précédents, et se convaincra que si la sonde effectue N tours complets avant le voyage classique à moins d'un tour, restant à faire, pour établir :

La résolution relèvera des mêmes calculs que plus haut.

2°) RESULTATS :

Cette fois la fonction F(a) n'est plus monotone , comme le montre l'exemple ci-dessous:

L'allure du graphe montre alors qu'une solution à 1 tour peut exister, à condition que le temps soit supérieur à un minimum, impossible à calculer littéralement ( du moins pour l'auteur ). Mais alors, au dessus de ce minimum, il existera 2 solutions à 1 tour.

Le propos pourrait se généraliser à plusieurs tours, sans grande difficulté.

VII RESOLUTION RAPIDE:

Le voyage est complètement précisé lorsque ayant caractérisé le cas de vol et résolu l'équation du temps qui donne a et b, on aura calculé les paramètres orbitaux de la conique. Donnons les résultats pour l'ellipse, laissant au lecteur le soin d'adapter le résultat au cas hyperbolique.

1°) Demi grand axe a et excentricité e :

C'est le calcul le plus facile.

Pour l'excentricité il suffit de partir des expressions de r₁ et r₂ en fonction de l'anomalie excentrique :

On tire alors aisément e :

Naturellement, l'excentricité e pourrait s'exprimer à l'aide de a et b ou de a uniquement.

3°) Vecteur vitesse :

La résolution du vol a permis de calculer pour chaque cas, a, b, Y, F, j₁, j₂, q₁, q₂, naturellement en apportant grand soin à la détermination correcte des angles.

NB : En particulier, l'inversion de j₁et j₂ avec les rayons vecteurs pose le problème de la détermination des angles, car la fonction y = Arc cosinus(x) donne un résultat 0 < y < p. Or il existe une autre solution entre p et 2p, y = 2p - Arc cosinus(x). Il faut donc vérifier laquelle des 2 solutions est la bonne, avant de continuer.

En utilisant les unitaires, radial et orthoradial, d'une position P1 ou P2, notés u et v, on vérifiera par exemple :

4°) Inclinaison orbitale i :

Le vecteur W qui oriente le plan orbital est parfaitement défini ( sauf cas particuliers et singuliers d = 0 ou p ) , et l'inclinaison orbitale i s'en déduit sans difficulté par :

5°) Longitude vernale W de la ligne des nœuds :

Sauf cas particulier i = 0° ou 180°, vraiment singuliers et pour lesquels classiquement W n'a plus de sens, un vecteur unitaire de la ligne des nœuds est bien défini par :

Donnant la longitude vernale de la ligne des nœuds :

6°) Argument nodal w du périgée :

Grâce au rayon vecteur r et à la vitesse V, on construit le vecteur excentricité e et tout en procédant à une vérification de e, on déduit l'argument nodal du périgée.

7°) Anomalie moyenne d'une date (départ ou arrivée) :

M = j - e*sinj calculée en P₁ ou P₂

7°) ROUTINES DE RESOLUTION :

Le problème de Lambert a été résolu informatiquement, pas nécessairement d'une manière qui convienne à tout le monde, mais utilisable. Deux langages ont été utilisés: Pascal et C++. Tous les cas sont envisagés y compris les vols hyperboliques.

a) La routine a été écrite en PASCAL par l'auteur du site et vous avez accès soit au source LAMBERT1.PAS, soit à l'exécutable LAMBERT1.EXE. Essayez pour voir.

b) La routine a été écrite en C++ par 2 étudiantes de la promotion DESS TAE 2002-2003 et vous avez accès soit aux éléments sources ( LAMBERTC.ZIP ) et au source LAMBERTC.EXE. Essayez pour voir.

RESOLUTION FINE DU PROBLEME DE LAMBERT

On conserve les conditions de base du problème de Lambert en faisant la différence entre une trajectoire qui relie 2 points en dehors de la notion de sphère d'influence d'une trajectoire plus réaliste, destinée à un transfert interplanétaire, ou encore un mélange des 2.

Ce qui caractérise le point de départ ou celui d'arrivée, c'est la sphère d'influence. Si elle a un rayon nul, c'est un point géométrique et le problème est parfaitement résolu. Si au contraire le rayon est non nul, il faut tenir compte de la phase de traversée de cette sphère d'influence.

I DONNEES

Nous notons R_spho et R_sph1 les rayons respectifs des sphères, avec possibilité de une ou deux valeurs nulles.

Nous conservons aussi D_o et D₁ les dates juliennes de départ et d'arrivée, prises au moment de l'évasion, c'est à dire à proximité de la planète, conventionnellement au périgée de l'hyperbole de départ ou de celle d'arrivée.

On suppose aussi que les éphémérides des planètes permettent de calculer Ro(t) =[ Xo(t), Yo(t), Zo(t)] et R1(t) =[ X1(t), Y1(t), Z1(t)] les coordonnées des planètes dans un repère héliocentrique bien précis, à tout instant.

II LES CONSEQUENCES DE L'OUBLI DES SPHERES D'INFLUENCE

La résolution du problème de Lambert est exacte entre 2 points géométriques.

Dans le cas de planètes, il apparaît que les vitesses à l'infini présentent de petites divergences notamment en ce qui concerne un tremplin, avec une évaluation de l'effet de tremplin faussée par la durée de traversée de la sphère, supposée instantanée, alors qu'elle ne l'est pas.

Illustrons par un calcul :

Le terme entre crochets est nul si on suppose un tremplin instantané, à t=t₁=t₂ , et on retrouve la formule classique du cours.

Par contre, la réalité indique que le tremplin dure un certain temps ( pour la Terre de l'ordre de 1 jour ), on peut dire que pour sortir de la sphère d'influence, très rapidement la sonde adopte sa vitesse de croisière, en s'éloignant sur l'asymptote à la vitesse V infinie( variant de 11.4 à 15.5 km/s). Le rayon de la sphère d'influence est sensiblement de 900000 km, ce qui donne un temps compris entre 0.67 et 0.9 jour.

Il est alors clair que par exemple pour la Terre avec une vitesse de 30 km/s, sur cette durée, le vecteur vitesse tourne de 0°.67 à 0°.9 environ, ce qui entraîne, prenons 1° une erreur maximale de 1 km/s :

L'erreur n'est pas tout à fait négligeable, entre 0.67 et 0.9 km/s.

III COMMENT Y REMEDIER

1°) AMELIORATION DE LA METHODE :

On calcule la trajectoire de Lambert joignant le point de sortie de la sphère d'influence de départ , à celui d'entrée dans la sphère d'influence à l'arrivée. Le départ est précédé d'une phase hyperbolique d'évasion, l'arrivée se poursuit par une phase hyperbolique de descente et de survol, avec ou sans tremplin gravitationnel.

Naturellement, la durée totale du voyage reste la même, mais celle du voyage de Lambert entre les sphères d'influence est plus courte et doit être évaluée finement. C'est l'objet des itérations au chapitre suivant.

2°) APPROXIMATIONS FAITES :

1 - Celle des sphères d'influence, dont l'expérience montre qu'elle est très bonne. Rayons notés Rd et Ra.

2- Les hyperboles de départ ou d'arrivée sont en pratiques tellement tendues que l'asymptote est à une distance très faible de 1 à 3 rayons planétaires du centre de la planète. Une telle distance permet de supposer que pour placer le point de sortie ou d'entrée dans la sphère, on peut le prendre de manière simple à partir du centre planétaire, sur la sphère, dans la direction de la vitesse à l'infini de sortie ou de son opposé dans le cas d'une rentrée. Un dessin montre mieux que tout cette configuration.

En clair les positions Po et Qo sont très voisines et sont donc confondues de même que P1 et Q1.

3°) METHODE NUMERIQUE :

On procède par itérations. L'indice de l'itération est l'indice du haut, celui du bas vaut d pour ce qui concerne le départ et a pour l'arrivée. On notera le vecteur position en limite de sphère d'influence, la vitesse héliocentrique en ce point de la sonde, celle de la planète, et la vitesse à l'infini pour la ième itération:

Initialisation: c'est le cas élémentaire ponctuel.

On résout le problème de Lambert et on trouve les vitesses héliocentriques, puis les vitesses à l'infini.

On calcule les temps de descente sur les hyperboles respectives de départ et d'arrivée, ces temps se nomment, avec des notations évidentes:

Rien n'empêche de calculer les temps exacts ( moins de 10% d'écart) , pour une précision certainement illusoire, mais peut être cela vaut la peine d'essayer. Il faudrait alors prendre en considération une hyperbole précise de descente, avec notamment une altitude de périgée imposée et conduire le calcul du temps de descente, avec l'anomalie excentrique hyperbolique.

Itération n : connaissant le niveau n-1

NB : Si on adopte le calcul simplifié du temps sur hyperbole, alors il vient :

Les vitesses à l'infini se calculent à chaque itération.

TEST D'ARRÊT DE L'ITERATION :

Vous devriez assister à une convergence de la solution, vous pouvez arrêter le processus par exemple en imposant :

ou peut être les 2 à la fois, c'est une question de pratique et de rapidité à tester.

GUIZIOU Robert décembre 2004, sept 2011

Existence d'un fichier Lambert.doc ( Word 97 ),

Avec une mise en page optimisée